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STABILIZED AND INEXACT ADAPTIVE METHODS FOR CAPTURING 
INTERNAL LAYERS IN QUASILINEAR PDE 


SARA POLLOCK 


ABSTRACT. A method is developed within an adaptive framework to solve quasilinear 
diffusion problems with internal and possibly boundary layers starting from a coarse 
mesh. The solution process is assumed to start on a mesh where the problem is badly 
resolved, and approximation properties of the exact problem and its corresponding fi¬ 
nite element solution do not hold. A sequence of stabilized and inexact partial solves 
allow the mesh to be refined to capture internal layers while an approximate solution 
is built eventually leading to an accurate approximation of both the problem and its 
solution. The innovations in the current work include a closed form definition for the 
numerical dissipation and inexact scaling parameters on each mesh refinement, as well 
as a convergence result for the residual of the discrete problem. Numerical experiments 
demonstrate the method on a range of problems featuring steep internal layers and high 
solution dependent frequencies of the diffusion coefficients. 


1. Introduction 


This investigation continues the work in f22l 211 developing adaptive numerical meth¬ 
ods for quasilinear partial differential equations featuring steep internal layers, in which 
the solution process starts on a coarse mesh where the problem is not yet resolved. 
Throughout the coarse mesh and preasymptotic regimes, standard methods such as New¬ 
ton iterations are known to fail due to both the ill-conditioned and possibly indefinite 
Jacobians which are characteristic of the approximate discrete problems, and the partial 
resolution of the problem data. This paper specifies an appropriate set of parameters that 
may be used in the stabilized cr-Newmark strategy of [1211 applied to quasilinear diffusion 
problems, where the layers develop from both the solution dependent coefficients and a 
variable dependent source. For the class of problems studied here 


— dw(tt(u) Vtf) — f(x ) = 0, u = 0 on and 

div(/s;(|V-u| 2 )Vw) — f(x ) = 0, u = 0 on dQ, 


( 1 . 1 ) 

( 1 . 2 ) 


in which the diffusion is bounded away from zero, local uniqueness of the solution is 
known, as well as approximation properties for the finite element solution using linear 
elements Il24l . assuming the mesh is sufficiently fine. For operators containing steep 
solution-dependent layers in the coefficients, the approximation properties of the discrete 
solution are useful only if the solution to the discrete nonlinear problem can be attained, 
and the current method attains such a solution by means of a sequence of approximate 
problems with inexact source functions that limit to the discrete problem. 

The methodology is to first discretize (If .f 1) and (11.21) on each mesh refinement then 
linearize the resulting discrete problem, and to partially solve the sequence of resulting 
inaccurate and ill-conditioned coarse mesh problems by stabilized Newton-like iterations 
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while adaptively refining the mesh leading to an accurate and efficient solve of an accu¬ 
rate discretization of the problem. In previous work the focus was on the stabilization of 
the Jacobian by a combination of regularization and added numerical dissipation. Cur¬ 
rently, a formula for the numerical dissipation parameter is presented, along with a new 
inexact method designed for problems where the variable-dependent source dominates 
the residual of the Newton-like iterations. 

Iterative rescaling techniques in the solution of nonlinear problems are not uncommon, 
see for instance the scaling iterative algorithm (SIA) of [J5]] for the solution of semilinear 
elliptic problems, in which the solution u is rescaled at each iteration. The rescaling 
of the Monge-Ampere equation in [|T| to establish a fixed-point argument and recover a 
numerical solution without having to assume the solution is small enough motivated the 
current approach. Here, the inexact method rescales the variable dependent source until 
the solution iterates attain sufficient stability to solve for the given data. 

Recent approaches such as flUlJT] f° r monotone quasilinear problems use inexact lin¬ 
ear and nonlinear solves to avoid over-solving for the residual when the Galerkin or 
discretization error are the dominant sources of error. It is assumed it their analysis that 
the discrete problem on each refinement is well posed. In the problems studied here, the 
coarse mesh problem may not be well posed, and may be a sufficiently bad approxima¬ 
tion of the exact problem that estimates of the different error sources are not necessarily 
well determined or useful. Instead, the iterations are ended when they stabilize to the 
predicted linear convergence rate, which is a function of the numerical dissipation pa¬ 
rameter. Combined with the criteria that the residual from the linear solves on each 
mesh refinement must show sufficient decrease with respect to the residual on the pre¬ 
vious refinement, the sequence of stabilized and inexact problems recovers the unsealed 
discrete problem on a mesh where it is better represented, and with an initial guess for 
the Newton-like iterations that yields the discrete problem solvable. This method pre¬ 
dicts the stability of the solve and allows the sequence of coarse mesh problems to be 
solved approximately through the preasymptotic regime leading to an efficient solve in 
the asymptotic regime. 

The remainder of the paper is organized as follows. Section [2] states the target prob¬ 
lem class and the formulation of the discrete problems. Section [3] reviews the Jacobian 
stabilization techniques developed by the author in previous work, and which are further 
developed here. Section |4]presents a formulation for the numerical dissipation parameter 
7 and characterizes its properties within the adaptive framework; then Section[5]presents 
a formulation for the inexact scaling parameter 5 and characterizes its convergence to 
unity within the adaptive method. Section [6] summarizes the results of the previous three 
sections into an adaptive algorithm and proves the convergence of the residual of the dis¬ 
crete problem. Finally, Section [7] demonstrates the method with a collection of numerical 
experiments featuring different types of internal layers. 

The following notation is used throughout the rest of the paper. In defining the weak 
and bilinear forms in the next section (u(x),v(x)) = f n u(x)v(x) dx, and in later section 
the discrete inner product between vectors Uk,Vk G M n is denoted (uk,Vk). The norm 
|| ■ || where not otherwise specified is the L 2 norm. The nth iterate subordinate to the A th 
partition Tk is denoted while u n is the nth iteration on a fixed partition and Uk is the 
final iteration on the fcth mesh, taken as the approximate solution on Tk- 

2. Target problem class 

The class of problems considered are quasilinear diffusion problems F(u, x) = 0, over 
polygonal domain Q C M 2 , with F : X xO —> Y* and F'(u, x) := F u (u, x) G C(X, Y*), 


LAYER CAPTURING ADAPTIVE METHODS 


3 


where F(u,x ) is given by 

F(u, x ) := —di y(k(u)'Vu) — f(x) — 0, in G M 2 , u = 0 on dfl, or (2.1) 

F(u, x) — div(/c(| Vw| 2 )Vtt) — f(x) =0, in Q G M 2 , u = 0 on <90, (2.2) 

with f(x) G L 2 (0) nL oc (fl). Multiplication against test functions u G F and integration 
by parts yields the weak form of each problem 

B(u,v) — (k(u)Vu, Vv) — (/, v), for all l G F, for (12.11) . (2.3) 

B(u, v) = («(|’ Vu\ 2 )Vu, Vv ) = (/, u) for all v G F, for Q). (2.4) 

The linearized form induced by F'(u, x) := F u (u, x), is determined by taking the Gateaux 
derivative in direction w G X by B'(u; w , v) = linp_^ 0 d/dt ( B(u + tw , u)) yielding 

B\u;w,v ) = (k(u)Ww, Vv) + Vr), for (12.11) . 

(2.5) 

B'(u]w,v ) = («(|Vm| 2 )Vr;, Vn) + (ac / (|Vm| 2 )(2Vm ■ Vtr)V«, Vn), for (12.21) . 

(2.6) 

Both types of problems fit into the context of [[24| with the assumption that there is a 
solutions G Hq(Q) D lF 2 2 + £ (^) an d F u (u,x) : Hq(Q) —* ff _1 (fl) is an isomorphism, 
in which case the solution u is an isolated solution, and approximation properties for the 
linear Lagrange finite element solution can be shown to hold, assuming the meshsize is 
fine enough. 

It is further shown in [@0] for problems of the form (12.11) . assuming n(s) > k > 0 for all 
sGl and «(s), k'(s), k"(s) bounded, then there exists a unique solution u G kL 1,p (fl), 
with 2 < p < oo. Moreover, by the analysis of [|T6ll . an adaptive method using linear 
Lagrange finite elements can be shown to converge, again assuming the meshsize is 
sufficiently small. 

The second class of equations (12.21) is considered in for instance [fl4l 3] |H and the 
reference therein, and uniqueness of solutions as well as convergence of adaptive meth¬ 
ods have been established supposing a Lipchitz condition and a monotonicity condi¬ 
tion such as B(v,v — w) — B(w,v — w) > 0 or the strong monotonicity condition 
B(v,v — w) — B(w,v — w) > c\\v — w\\y as in ItTO . and assuming the approximate 
discrete problems are well posed. 

2.1. Discrete formulation. The discrete problem is, find Uh G Xj, such that /i (»/,., v) = 
0 for all v G Y h where X h C X and Y h C Y are discrete spaces of the same finite 
dimension. Here, X h and Y h are finite element spaces with respect to triangulation 
where the family of triangulations {Th}o<h<i is regular and quasi-uniform in the sense 
of 0. For the rest of the discussion it is assumed the trial and test spaces are the same: 
X = Y and Xj, = 1),. On each refinement k, the test/trial space Xh k is taken as the 
linear finite element space 14 consisting of Lagrange finite elements Pi over partition % 
that satisfy the homogeneous Dirichlet boundary conditions. 

The resulting discrete problem on refinement k is of the form Gk(u- x) — 0 for u G 14 
and x G H C M 2 . The decomposition Gk{u,x ) = — g(u ) + fk(x) is used throughout 
the remainder of the discussion. On a given discretization that is not sufficiently fine, 
fk(x) may be a bad representation of the given source f(x). On such a coarse mesh, 
solving the discretized equation g{u k ) = f k {x) may produce a solution u k which is a 
highly inaccurate solution to (12.31) . or respectively (12.41) . so seeking an exact solution to 
the approximate problem is not a good use of computational effort. 
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Ideally, the adaptive method starts on a coarse mesh where the problem is small, if 
badly resolved. It produces a sequence of approximations that are good enough to refine 
the mesh to resolve the data and solution dependent coefficients, and upon sufficient 
resolution and stability, provide a good approximation for the solution to the discrete 
problem. 

The features of the current method as introduced in If22l t21l include an adaptive ap¬ 
proach to regularization used to stabili z e the Jacobian in the coarse mesh and preasymp- 
totic regimes, also conceived of as adaptively regulated implicit pseudo-time stepping. 

3. Stabilized adaptive method 

The stabilized iterations are designed to fit into the framework of a standard adaptive 
method, which uses error indicators to mark and refine the mesh. The theoretical results 
of this paper are independent of the specifics of the adaptive method, and a different set 
of error indicators could be used. The method here is demonstrated using standard resid¬ 
ual based estimators, as in [12] [23] for mesh refinement, and similarly for the adaptive 
regularization strategy. Both the interior and the jump part of the residual indicator are 
used to refine the mesh, and a similar jump indicator based only on the solution is used 
to determine the regularization. The indicators are defined with respect to the consistent 
problem with unsealed data. 

3.1. Adaptive mesh refinement. The local a posteriori residual-based indicator, and 
corresponding jump-based indicator for element T E Tk with /; 7 the element diameter 
are given by 

Ct( v ) = CnM '■= hr\\ Jriy') IIl 2 (9t)j (3-1) 

Vt(v) = Vr k ( v ,T) '■= h T\\ F ( v , x )\\l 2 (T) + (t( v ), ( 3 - 2 ) 

Jt{v) ■— |[k(p)Vp • njor, for problem (12.11) and J r (n) := [[k(|Vp| 2 )Vp • njgr, for 
problem (12.21) . with jump : = lim t _> 0 4>(x + tn ) — f(x — tn ), where n is the appro¬ 
priate outward normal defined on dT. The error estimator on partition Tk is given by the 
l 2 sum of indicators rjj- k = ^2 T&7 - k ifr, and similarly for Qr k ■ 

In these results, the Dorfler marking strategy is used. Given a parameter 0 E (0,1) a 
set of least, or nearly least cardinality is chosen for the marked set M. so the sum of their 
indicators is greater than the given fraction of the sum of all the indicators. 

^Vt> 0 Y1 r &- ( 3 -3) 

TeM TeTfc 

Remark 3.1. In contrast to an earlier version of the method for convection diffusion 
problems in fl2Tll . no additional coarse mesh marking is performed. Scaling down the 
variable dependent source on the initial solves cdlows the method to traverse the preasymp- 
totic and coarse mesh regimes while minimizing refinement that negatively impacts the 
efficiency of the method in the asymptotic regime. 

3.2. Adaptive regularization. The motivation from multiple perspectives behind the 
adaptive regularization technique from is described in lf22ll . and an updated strategy is 
suggested in ll2Tll . which is summarized below. The coarse mesh problems may have 
indefinite, and often have highly ill-conditioned Jacobians, and the Laplacian-based reg¬ 
ularization helps make the Newton-like iterations solvable. The adaptive element is in¬ 
troduced to yield faster convergence as the degrees of freedom for which the approximate 
solution is smooth experience little benefit from the added stabilization as compared to 
the degrees of freedom (dofs) in regions where the jump in the normal derivative between 
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neighboring elements is large. On each mesh partition 7/ the penalty matrix R = R k is a 
localized version of the Laplacian stiffness matrix denoted Ra lobal _ First define the jump 
indicator used to select elements for regularization, c.f, (13.11) . 

€t( v ) = M |[Vu • nj dT ||i a(ar) . (3.4) 

Using the jump indicators select degrees of freedom for regularization. 


Definition 3.2. Define R' llohal as the Laplacian stiffness matrix on refinement k. Let 
D loc a diagonal matrix of zeros and ones, v 3 a vertex subordinate to partition Tk, and 
Then set 

fk = \f• median T( zj- (ff ), and 'iff — { ^~ k ’ > 

V V | ^ otherwise, 


and 


D 


loc 

33 


0 , if ff < V'k for each element that contains v 3 as a vertex , 

1 , otherwise. 


(3.5) 


Set R k = D loc RP lobal D loc . 


The regularization matrix R k defined by Definition 13.21 is positive semidefinite and 
targeted to smooth regions of high curvature. 


3.3. Exact and inexact cr-Newmark iterations. The cr-Newark iteration of [[2TTl gen¬ 
eralizes the pseudo-transient continuation method of [18] [ 8 , 13] and similarly the “s” 
method of [[ 2 ]| to an adaptive framework, combining stabilization by pseudo-time with 
the Tikhonov-type regularization as used in ill-posed problems ITOTl . 

The Jacobian is stabilized by seeking a steady-state solution of the given nonlinear 
problem Gffupx) = 0 by solving d(Ru)/dt + G k (u,x ) = 0, and discretizing the time 
derivative u := du/dt by a Newmark update |[20|| given by the time discretization v n+ 1 = 
u n + A t n { n fu n+1 + (1 — 7 )u n }, or solving for ii n+1 


u n+1 -u n ( 1 - 7 ) , 

--- u 

7 A t n 


In the current notation the regularization parameter a n = 1/A t n . As discussed in the 
literature on the finite element analysis of structures, for instance [1151 [ 6 l and the refer¬ 
ences therein, the update yields second-order accuracy in time for 7 = 1/2 and increased 
numerical dissipation of the high frequency modes with 7 > 1/2. The current adaptive 
method allows 7 > 1 , which can be seen to damp the influence of the right-hand side 
residual. In [[211 a secondary linearization is performed to further stabilize the approxi¬ 
mate Jacobian resulting in the iteration 


(a n R + 7 n {a n G’{u n ) + (1 - o n )G\u)}) w n 


—G(u n , x ), u n+1 = u n + w n . (3.6) 


In previous discussions the variable dependence of G(u, x) was suppressed and written 
G{u). In the current presentation the variable dependent load will be stated explicitly 
by G(u,x) = g(u ) — f(x). As the Jacobian G u (u,x ) = g’(u ) is invariant to f(x), the 
stabilization based on the Jacobian alone does not address sharp spikes or steep gradient 
which may be present in f(x), particularly in partially resolved coarse mesh problems. 
The inexact iteration is then introduced 


(, a n R + 7 n {o n g\u n ) + (1 - affgfu)}) w n = -g(u n ) + 5f(x), u n+1 = u n + w n . 

(3.7) 
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Iteration (13.61) is defined by the following parameters on iteration n of refinement k. The 
limiting behaviors necessary to reduce the stabilized iteration to the standard Newton 
iteration g'{u)w n = — g(u n ) + f(x) are noted. 

• Regularization parameter a = a k > 0, a — Y 0. 

• Semidefinite regularization term R = R k . 

• Newmark parameter 7 = 7^ > 1, 7^ —> 1. 

• Stabilization parameter a = a k G [ 07 , 1], a —> 1. 

In addition, iteration (13.71) depends on the scaling parameter 

S — 8k < 1, 5 —» 1. 

The focus of this paper is the description of an appropriate choice of parameters 7 and 5 
to use in iterations (13.61) and (13.71) . The regularization parameters a and R are presented 
in Itm and recalled in this section. The convergence behavior of iteration (13.61) is ana¬ 
lyzed in [f2Tl . and those results are recalled below. Definitions of 7 and S are introduced 
in the following sections, as well as a convergence result for the sequence of inexact 
iterations using (13.71) . 

The following concepts are used throughout the discussion. On a fixed refinement 
k, the subscript k is removed where it will not cause confusion, e.g., the discretized 
variable-dependent source f k is denoted /. The solution dependent part of the residual 
g(u n ) is denoted g n , and residual r n is given by 

r n — ~g n + Sf, (3.8) 

where S = 1 for iteration (13.61) . 

Before introducing the definitions of the numerical dissipation parameter 7 and the 
scaling parameter 5, local convergence the a-Newmark iteration (13.61) is recalled from [fTTI 
As the scaling parameter <5 is held constant over iterations on a given refinement, The¬ 
orem 13.61 holds also for the inexact iteration (13.71) . where in the first case the residual 
which converges to zero is given by r(u, x ) = —g(u) + f(x) and in the second it is given 
by r(u, x) = — g{u) + Sf(x). The proof of local convergence depends on the nonlinear 
discrete problem, which coincides up to sign with the residual function on a fixed re¬ 
finement, satisfying a local Lipschitz condition on the Jacobian with respect to u in the 
vicinity of the exact solution u*. For simplicity of notation, consider the convergence of 
—r(u, x) to zero, whose Jacobian is then g'{u). Denote the open neighborhood about u 
by N(u, e) = {v | ||m — v\\ < s}. 

Assumption 3.3. There exist u L ,e > 0 so that for all w,v G N (u*\ e) 

|| g'(w) — g\v)\\ < oolW'w ~ f|| for all w, v G N(u*, e). 

The second set of assumptions addresses the nonsingularity, boundedness and stability 
of the approximate Jacobian. Assumption 13 .51 requires a regularization parameter a n < 
OiM for some fixed om, which follows from the definition of o” described in ll22i and 
repeated here for convenience. In accordance with the local convergence Theorem 13.61 
the update assures a n < ||r n ||, so long as the residual is decreasing. 

Definition 3.4 (Regularization parameter a). Set /3° = 1. For n > 1, 

a n = /3 n ||r n ||, with 

on _ \ min{l ,max{/3 n ~ 1 /2 , ||L n ||/||r n-1 ||}}, (/’||r n || < ||r ?l_1 ||, 

1 f> n ~ l 1 otherwise. 
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The factor /3 n of Definition 13.41 modifies a n to reduce the regularization faster if the 
residual is decreasing steadily, without inducing rapid fluctuations in the regularization 
between iterations. So long as the residual is decreasing, a n < ||r n || should be main¬ 
tained to satisfy the hypotheses of the local convergence Theorem 13.61 The last set of 
assumptions for the local convergence theorem addresses the invertibility and stability of 
the approximate Jacobian. 


Assumption 3.5. There exist 3 > 0 so that for 0 < a 0 < o < 1, fixed u, positive 
semidefinite R, 7 > 1’ and for all 0 < a n < <Xm, then for all u G N(u*,e): 

1) a n R + 7{(1 — o)g'{u) — crg'{u)} is invertible. 

2) \\{a n R + 7{(1 — o)g'{u) — cr(/(w)}) _1 || < M 1<7 . 

3) \\(a n R)(a n R + ^{1 - a)g'(u) - agfu)})- 1 ]] < 1+p \ /an - 

The local convergence theorem demonstrates that for a sufficiently small residual, if 
the incoming iterate u n is within the basin of attraction of the exact solution ?/*, then u n 
is also within the basin of attraction, and the residual decreases asymptotically at the rate 

|| r n+i|| < ^ _ l/y)||r ri ||, for 7 > 1. 


Theorem 3.6. Let a n < ||r(w ri )|| < an d let Assumptions [PI and \T5\ hold. and a n 
given by Defniteion \3 A\ Define a by 


o = max < cr 0 , 1 


/In 


(3.9) 


for a given K 0 . Then there exists d > 0 such that for v ri in the open set S given by 

5 = («e N(u*,e) ||r(w)|| < d\ , (3.10) 


and u" + 1 defined by iteration (13.61) . or respectively by (13.71) with 7 > 1, it holds that 
G S, and the sequence of residuals converges q-linearly to zero with asymptotic 


u 


n +1 


rate q = 1 — I/7. 

The proof of Theorem I3T61 is given in [12TI . 


4. Update of numerical dissipation. 

The numerical dissipation parameter 7 > 1 of iterations (13.61) and (13.71) should start 
large enough to provide sufficient stability to the iterations and decrease down to one 
in order to recover the quadratic convergence rate of the Newon-like iterations. The 
parameter is adjusted during the course of iterations on a given refinement whenever the 
observed rate of convergence is within tolerance of the expected rate, (1 — I/7). In order 
to ensure a decrease in 7 as it is updated, each 7]? must be bounded by an appropriate 
7 m AX , with respect to the user set rate tolerance, e T - 

Assumption 4.1. Given a rate tolerance Et, LS bounded above for all n, k by 7 max > 
satisfying 

7 k < Imax < 1 /e T . (4.1) 

This assumption has been observed to be necessary, and 7 can indeed grow without 
bound if allowed greater than 1 /e t . It is also of practical value not to set the rate tolerance 
Et too low, as the update of 7 can stall for some number of refinements due to the partial 
capturing of the solution and variable dependent layers in the preasymptotic regime. In 
the examples presented in Section |7j e t is chosen either as 0.005 or 0.01, where the 
larger value is chosen in the case where steep gradients prevent the resolution of the 
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Jacobian, and this adjustment is sufficient to speed the update of 7 without compromising 
stability of the method. With these parameters, the update of 7 on a given refinement k 
is performed when the following criteria are met. 


Criteria 4.2 (7 update). Given a user set tolerance c T , and 7"’ < 7 max satisfying As¬ 
sumption [O] and a minimum number of iterations between updates I m j n , the dissipation 
parameter y n is updated after each iteration n on which all the following criteria hold. 


1 ) 7" > 1. 

2) The observed rate of convergence is within tolerance of the predicted value 


„n+l I 


- 1- 


ry T > 


< £t- 


3) The observed rate of convergence is sufficiently stable. 


„n+11 


r.n—1 


< Ef. 


(4.2) 


(4.3) 


4) At least I min > 2 iterations have passed since the last update of y. 


The last condition of Criteria l4~2l is necessary to assure the three consecutive residuals 
of (14.31) are computed with the same value of 7, and show the requisite rate stability. 
Given a reduction factor q — q 1 < 1, on satisfaction of Criteria l4~2l after iteration n of 
fixed refinement k, y n+1 is updated as follows. 


Definition 4.3. Given an initial 5 0 > l, for k > 0 

(r n ,r n ) 


y n+1 = q ■ 


(r n , g n+1 - g n ) 


, y n+l = max {1 , y n+1 } . 


(4.4) 


To see this update is reasonable as the algorithm attains convergence, observe the ideal 
value of g n+l is Sf, implying r n+1 = 0, and the ideal value of g n+1 — g 11 is Sf — g n = r n , 
so that y n+l = 1, if the iteration produces a zero residual. Based on condition (14.21) of 
Assumption 14.21 the updated quantity y n+1 is bounded above and below in terms of y n , 
and is then shown to be decreasing, on the condition of Assumption 14.61 Lemma l4~4l 
bounds the denominator (r n , g n+1 — g n ), then Corollary 14 .5 I bounds the ratio determining 

7 n+1 . 


Lemma 4.4. Under (14.2b of Assumption 1 4~2\ with tolerance £t the following upper and 
lower bounds hold. 


- n " 2 1 1 -£ T ) < (r n ,g n+1 - g n ) < \\r n \\ 2 ( 2- 1 


ryi 




+ £t 


Proof. Rewriting g n and g n+l in terms of the residual 

(r 

For the first inequality 


r“, j” +1 - 9 n ) = (r“, -g" + Sf- (-j” +1 + Sf)) = (r”,r' 


n+1 


)• 


(jn _ (y.n r n +^ ||r n || 2 — ||r"|||| r n+1 1| = ||r n || 2 ( 1 


„n+11 


(4.5) 


(4.6) 


(4.7) 


Bounding the ratio of residuals on the right by (14.21) 

llr. n+l II 


1 - 


> - — £t. 

ryTL 


Applying (14.71) and (14.81) to (14.61) . the result follows. 


(4.8) 
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For the second inequality of (14.51) . applying Cauchy-Schwarz and triangle inequalities 
to (l4~6l) 

( j> n r n _ r n + 1 '^ < lumi || r ” _ r n+1 || < || r n || 2 _|_ ||r n || ||r n+1 1 | (4 9 ) 

By (14.21) . ||r n+1 || < ||r n ||(l — l/y + £t), from which the result follows. □ 

The upper and lower bounds on 7 n+1 now follow from (14.51) . 


Corollary 4.5. Under Assumption ^. h and Criteria \4~2\ the update of y n to y" 4 1 of (14.41) 
satisfies the following bounds. 


q . 

y n (2 + e T ) - 1 


< 7 


n+1 


< 


q ■ y n 
1 — CtT™ 


(4.10) 


Proof The argument follows from direct application of Lemma 14.41 to the definition of 
7 n+1 given by (14.41) . The upper bound 

11^,71112 1 


7 n+1 = q • 


< q- 


q * 7 


{r' b , g n+1 — g n ) 1 / 7 n — 1 — £r 7 n ’ 

holds under Assumption 14. 11 which assures the positivity of the denominator. Similarly, 
the lower bound follows applying 7 > 1, the first condition of update Criteria 1431 □ 


The lower bound remains satisfied in the case 7 n+1 = 1 by Definition 14.31 and the 
upper bound is of no interest as the goal of the sequence is to decrease 7 to unity in a 
stable manner. To establish the monotonic decrease of 7, and more generally the decrease 
at rate 1 > q > q, the following Assumption l4.6l strengthens the Assumption l4.ll used to 
establish the upper bound of (14.101) . 


Assumption 4.6 (Monotone decrease of 7). Given a reduction factor q, a rate of decrease 
q > q, and rate tolerance £t, let 7 ^ be bounded above for all n, k by 7 mono satisfying 


7 k — Imono 



(4.11) 


In practice, setting q — 1 and 7 < (1 — cf/ej ensures 7 will be nonincreasing, although 
for more stability in the earlier iterations, the weaker condition 7 < 7 max is sufficient, 
and after 7 reduces below 7 mono the monotonicity holds. 


Corollary 4.7 (Monotonic decrease of 7). For y n+1 determined by (14.41) on the satisfac¬ 
tion ofCriteria \4.2\ and Assumption \4 .61 with q — 1, the sequence {y n }n>o is nonincreas¬ 
ing. Moreover, assuming (14.1 II) with 1 < q < q, then for y n+1 = y n+1 > 1, the sequence 
satisfies y n+l < qy n . 


The proof follows from the upper bound on 7 n+1 , given by (14.101) of Corollary 14.51 

Proof. To assure 7 n+1 < qy n /(I — £ry n ) < qy n , require 

q < q(l- £ T y n ), (4.12) 

which is satisfied by Assumption 14.61 To establish the monotonic decrease without a 
rate, set q = 1 in (14.121) . and rearrange terms to obtain 7 " < (1 — q) je-p as the necessary 
condition. □ 


While the determination of the monotonic decrease of 7 with a given rate q to unity 
follows in a straightforward manner from Definition 14.31 condition (14.21) . and Assump¬ 
tion [431 this establishes the important property that if the algorithm does not reset after 
a given point, 7^ = 1 will be achieved after at most N — n= \ — log(7 n )/ log(g)] up¬ 
dates, where [•] denotes the ceiling function. Notably, N — n counts only the updates 
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of 7 , independent of the number of mesh refinements that occur over the course of those 
updates. The more delicate point of whether those updates occur, which requires the pre¬ 
dicted convergence rate is achieved within tolerance, is addressed in Theorem 16.61 which 
relies on Theorem 13 .61 the local convergence with rate of iterations (13.61) and (13.71) . 


5. Update of scaling parameter 


The scaling parameter <5 for the inexact iteration (13.71) is computed on the completion 
of the iterations for refinement k. This is in contrast to the update on 7 , performed 
whenever the iterations converge close enough to the predicted rate. This strategy is 
reasonable as updating 5 changes the discrete problem which already changes on each 
refinement. Except for the first, the criteria to update 5 are a subset of the conditions 
used to exit the Newton-like iterations on each refinement. The last two conditions of 
Criteria 15.11 are closely related to (14.21) and (14.31) . the second and third conditions of 
Criteria 14.21 the update on 7 . The acceptable rate criterium (15.21) is weaker than the 
requirement (14.21) . and the stability criterium (15.31) requires the observed rate has not 
decreased by much compared to the previous observed rate. Distinct from the conditions 
to update 7 , the second condition of Criteria l5Tl requires the current and previous residual 
to be smaller than the initial residual of the current refinement and final residual of the 
previous refinement. This condition is important to ensure the algorithm convergences, 
and the sequence of residuals r k decrease at a rate 

INI < lkfc-ill(i - 1 /( 27 *)) < Ikfc-ilK 1 - V( 2 7max)), 

on Assumption l4.ll or 

INI < ll r fc-i||(l - l/( 2 7fc)) < IlUc—i||(l - 1/( 2 7 mono)), 

on Assumption 14.61 where 77 is the final residual on refinement k. Conditions (15.21) 
and (15.31) are also important to the stability of the update, and assure g n+1 = g{u n + w n ) 
has been improved by two consecutive steps w n i and w n . 


Criteria 5.1 (8 update). Given a user set rate tolerance Et as in Criteria 14.21 and a 
convergence tolerance Econ, the iterations are terminated and 8 is updated after iteration 
n of refinement k, when conditions (1-4) are all met, or when conditions (1) and (5) are 
satisfied. 


1) 8 k < 1 . If 8k = 1, it is not updated. 

2) Sufficient decrease in the sequence of residuals 

Ikfcll < min (lNl, INi||}. 

3) The obserx’ed rate of convergence is beneath the accepted rate (1 — I/ 27 ) 

ll^.n+ 1 1| 1 


< 1 


27^ 


4) The observed rate is sufficiently stable 

||,v,nTl || 


+ ^> 


r n—l | 


(5.1) 


(5.2) 


(5.3) 


||T* ,£, || 2 W^' 

5) 11 r” | j < econ■ The iteration reached convergence. 

On satisfaction of Criteria 15711 8 is updated based on replacing the satisfied linearized 
equation 

{a n Rk + l n (o- n g\u n ) + (1 - cr n ) g '(u))}w n = -g n + 8 k f k , (5.4) 
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with the corresponding relation found by replacing g'(u n )w n , with g n+l — g n , modifying 
the right hand side residual by q < 1 to obtain 

{a n R k + 7 n ((1 - a n )g'(u))}w n + 7 V" ( g n+1 - g n ) = -g n + qS k+1 f k , (5.5) 

and considering the discrete inner product of both sides against the source / = f k (x). 
Given a factor q = qs < 1, at the end of each nonlinear solve, 4+i is determined by 

{/, {c n R t + 7 " ((1 - + 7 V" ( 9 " +1 - j”)> = (/, —g n + q t S t+1 f), 

(5.6) 

resulting in the following definition. 


Definition 5.2. Given an initial 4 < 1 , for k > 0 

(/, 7 °{g n+1 - g n ) + ( 7(1 - cr)g'(u) + aR)w n + g n ) 


4 +1 — 


4+1 = min{l , 4 +i}, 


Qk\\f\\' 


(5.7) 

(5.8) 


with R = R k , a = a n , 7 = 7 ” cr = <j n . 


The sequence for 4 < 1 can also be expressed recursively by considering (15.41) . satis¬ 
fied exactly with 8 = 8 k , and subtracting the corresponding inner product against / from 
(15.61) . Expressed in terms of the one-step linearization error C e {u n ) given by 

£>”) := S ” +1 - 9 “ - </(«’>”, 9 ” +1 = 9 («” + w n ), (5.9) 

the recursive relation for 5 may be written 

4 +1 = 4 E + ^j_ (/i £=(„»))|. (5.10) 

From (15.101) . the local Lipschitz Assumption 13 .3 1 on g'(u k ) is used to show the conver¬ 
gence of 8 to one as the residual attains convergence. First, a bound on 4+i with respect 
to 4 is established. 


Lemma 5.3. Assume u k G N{u* k: e k ) of the local Lipschitz condition Assumption 13.31 
with respect to the approximate Jacobian g'(u k ), and where u* k satisfies r k {u k ) = 0. 
Let uiL,k denote the corresponding Lipschitz constant. Then 4+i updated by (15.71) with 
q k < 1 , satisfies the following bound with respect to 8 k . 


1 

qk 



T aU} L,k 

2H/II 



< 5 



Qk 



7 

2H/II 



(5.11) 


Proof. Writing 4+i by the recursive definition given by (15.101) . writing out C e (u n ) ex¬ 
plicitly 

4+i = 4(4+ g(u n+1 ) - g(u n ) - g'KK*}} , (5.12) 

Qk l \\f\\ 2 ] 

with u n+1 = u n + w n , yielding from the integral mean value theorem 

g{u n+l ) - g(u n ) - g'(u n )w n = [ {g\u n + tw n ) - g\u n )) w n dt. (5.13) 

Jo 

Applying a Cauchy-Schwarz inequality to (15.131) . and the Fipschitz Assumption l3.3l 

(fj\g\u n + tw n )~g'(u n ))w n dt\ < \\f\Pf\\w n \\ 2 . (5.14) 
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Putting together (15.131) and (15.141) into (15.121) yields the first result, (15.111) . This shows 
where the Lipschitz estimate on Jacobian holds, the variance in <5 fc+1 from a strict reduc¬ 
tion by factor q k is bounded with respect to the Lipschitz constant uiL,k- □ 

The next Corollary 15.61 establishes the monotonicity of the sequence {S k }k>K 0 for 
some finite K 0 , assuming the sequence of residuals ||r fc || is decreasing. This is imposed 
by assuming 5 is updated on the satisfaction of Criteria 15711 meaning the iterations con¬ 
verge within the accepted rate on each refinement. An easily imposed lower bound for 5 
is assumed. 


Assumption 5.4. There is a constant 0 < 8 M in < 1 for which 0 < 8 M in < A: for all 
k > 0. 


This assumption can be imposed in the definition of 8 k+ i from 8 k but is given sepa¬ 
rately for clarity of presentation. The lower bound <5 MIN is related to generally inaccessi¬ 
ble constants which characterize the approximate discrete problems in the next Assump¬ 
tion [531 As a guideline, 3 M in = 1/7max is reasonable as 8 tends to increase roughly 
as l/y. A need to impose 8 > 3 M in was not encountered in the numerical experiments 
presented in Section [7] Along with 15 .41 a condition on the data of the approximate prob¬ 
lems is assumed, requiring the ratio of 0JL,k, the Lipschitz constant for the approximate 
problem on refinement k, to \\fk\\, remains bounded. 


Assumption 5.5. There is a constant B^f for which after some refinement K\ 

for all k>K u (5.15) 

ll/fell 

and the convergence tolerance e C on of Criteria 15. l\ satisfies 


2 ^ 8min{ 1 Qs/q) 2 1 , 

£ con < -33;- T7f° raU k > A i- ( 5 - 16 ) 

Ik a k 1V1 ~ia 

Corollary 5.6 (Convergence of 8 to unity). Let Assumption 13.51 and the hypotheses of 
Lemma 15.31 hold. Then b k+ \ updated on satisfaction of Criteria 15.71 with q k < q$ < 1 
satisfies 8 k+ 1 > 8 k /q for q$ < q < 1 whenever 


7 V ra Ikl 2 M 2 ia u L ,k 

(1 ~qs/q) 2 H/fcll 


(5.17) 


Additionally, if Assumptions 15.41 and 15.51 hold, then (15.171) holds for all refinements k > 
K 0 > K x , so long as the residual ||r^ 0 || satisfies 


I K II 2 < ^ MIN if — Qs/q) 


,-w 77.0-71 


M 2 

7cr 


B zh 


(5.18) 


and 8 k = 1 for all k > K for some K < flog(A jFi - 0 )/ logg] + K 0 , where [■] denotes the 
ceiling function. 


Proof. Denoting the approximate Jacobian on iteration n by 


M = ( a n R + 7 n ( ag'(u n ) + (1 - a)g'(u))). 

In accordance with Assumption 13.51 ||M -1 || < M ia , and ||w n || < ||r n ||M 7(T . Then 
applying q s > q k 

T n cr n ||f n || 2 M> i|t N 7 n ° n \\w n \\ 2 u} L , k 

(1 -qs/q) 2 \\f k \\ ~ (l-q k /q) 2 \\f k W ^ J 
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Satisfaction of (15.171) then implies 5k > 5/q when dp. is greater than the quantity on the 
right of (15.191) . Rearranging terms in (15.191) and comparison with (15.1 II) yields the first 
result 

X 1 f x ^ aUL ^\\„..rm2 \ ^ 

Ofc+i > — 5k - I, ... w I > -bk- (5.20) 

Qk V 211/11 J Q 

Similarly for the second result, B~j < \\fk\\^L,k, 5min < 5k, and for k > K 0 , Crite¬ 
ria [5J] conditions either (2) or (4) and Assumption 15.51 yield \\rk\\ < ||rK 0 ||, by which 
satisfaction of (15.181) implies 


2 < 5k(l-qs/q) 2 \\fk\\ 

- y n a n to K)L ' 


(5.21) 


Rearranging terms again shows (15.171) . So far, this shows A /,. + 1 > 5k/q for k > K 0 . 
A short calculation and Definition 15.21 yield the final result, c) fc = 1 after at most an 
additional |~log(5x 0 )/ log <?] refinements and updates of 5. □ 


Importantly, the recursive definition (15.101) shows 5 is increased more aggressively 
when the linearization error lies in the direction of /, and is increased less and possibly 
decreased when C e (u n ) is in opposition to /. This sensitivity to direction as opposed 
to magnitude alone yields added stability to the method. More precisely, the following 
proposition characterizes the conditions supporting increase of 5 in terms of the projec¬ 
tion of the linearization error along /. The normalized vector in the direction of / is 
denoted / = ( 1 /||/||)/. 


Proposition 5.7 (Characterization of increase in 5). The sequence {^} given by the 
recursive definition (15.101) satisfies 5k +1 > (1 /q)5k, 1 < q < (p- whenever 

{/, £■(«")) > -4 (l - 4) (5.22) 

v 9 j 

The criteria (15.221) follows by direct calculation from (15.101) . The lower bound on (15.1 II) 
ensures the increase in 5k in the asymptotic regime as 11 w" | —y 0 on a given refinement 
assuring the algorithm eventually convergences to the given unsealed data, 5 = 1; how¬ 
ever, in practice the increase in is governed by (15.221) through the preasymptotic and 
coarse mesh regimes. 


6. Algorithm 

The results of the previous sections are summarized in the following algorithm, and 
combined to yield a proof of convergence using iteration (13.71) . The solver using inexact 
iteration (13.71) is implemented in an adaptive method according to the following basic 
algorithm. The algorithm using (13.61) is the same, except 5 = 1 is assumed throughout. 

The additional criteria and parameter settings necessary to define the method are sum¬ 
marized below the algorithm. The choice of parameters to guarantee the monotonic de¬ 
crease of 7 as in Corollary 14 .7 1 is listed as optional. The selection of parameters 5 M in and 
£con t° guarantee the monotonic increase in 5 as in Corollary 15.61 depend on unknown 
constants that are not assumed available. 

Algorithm 6.1 (Algorithm using the inexact iteration (13.71) ). Set the parameters g 7 , e T 
and 7 max in accordance with Assumption I4.il and optionally 14.61 Set the convergence 
tolerance £ C on and 5 M in of Assumption 15.41 Let u — 0. Start with initial u°, 7 0 , 5 (l . On 
partition %, k = 0 , 1 , 2 ,... 

1) Compute Rk by Definition ^. 2\ and g'[u). 
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2) Set r° = —g(u°) + 8 0 f 0 , and set a 0 = ||r°||. 

3) While Exit criteria \6.2\ are not met on iteration n — 1 : 

(i) Set a according to (13.91) . 

(ii) Solve ( a n Rk + 7 n {o n g'(u n ) + (1 — a)g'(u)})w n = r n , for w n . 

(iii) Update u n+1 = u n + w n , and r n+1 = -g{u n+r ) + S k f k . 

(iv) If Criteria \4 .21 are satisfied, update y n+l according to Definition \4.3\ 

(v) Update o" according to (13.41) . 

4) If Criteria 15. h are satisfied, update 8 k+ \for partition T k +i according to (15.21) with 
qs = (q-, ) p+1 , where P counts the number of times 7 was updated on refinement k. 

5) Compute the error indicators to determine the next mesh refinement. 


Criteria 6.2 (Exit criteria). Given a user set tolerance £con> an d an accepted rate of 
convergence given by q AC c = 1 — 1/( 27 ) as in Criteria I5.il and a baseline number of 
allowed iterations Ibase> set the maximum number of iterations Imax by 


Iacc — 


log(lk-i||) - log M|)' 


+ 1, Imax — max{/ ACC , I base}- 


( 6 . 1 ) 


log (q A cc) 

Exit the solver on partition T k after calculating residual ||r n+l | when one of the following 
conditions holds. 


1) Condition (5) of Criteria \5.1\ the iteration reached convergence. 

2) Conditions (2-4) of Criteria \5.1\ the iteration converges an at acceptable rate and 
is stopped before convergence. 

3) Maximum number of iterations I M ax completed: reset y k+l , 8 k+ \, andu k+l . 

Remark 6.3. Setting the maximum number of iterations bv \6.1\ prevents a premature exit 
when the residual is converging close to the predicted rate. An additional failure criteria 
may be added to exit the iterations quickly if the convergence rate drifts substantially 
higher that the accepted rate. However, this is mainly useful on very coarse meshes 
where the problem is small and the linear solves take negligible time. 

Remark 6.4 (Initial and reset values for 7, 3, u). An initial set of values for 7, 8, and u 
are not formally a feature of Algorithm 16. /I however some set of values is required to 
start the algorithm. The examples shown in Section\7\all use 

<=0 ' 7S= tik' ' o=1/tS ' <6 - 2) 

If the iterations fail to converge ,i.c., the iterations are terminated on iteration k by 
condition (3) of Criteria 16.21 u is reset and 7 and 5 are updated from their previous 
values by 

u ° k+ 1 = 0, 7° +1 = min{ 27 fc, 7 M 4 x}, 8 k+ i = max{8 k /2, 5 M in}- (6.3) 

Alternatively, 7 mono can take the place of "(max in (16.31) . 

Allowing the possibility of reset is considered an essential part of the algorithm, as 
the method is designed for capturing steep layers starting from a coarse mesh where the 
layers may be only partially uncovered or even completely undetected. The approximate 
discrete problem can drastically change over certain mesh refinements, and the latest 
update of parameters based on a milder approximate problem may not provide sufficient 
stability. 

The next Assumption[63] supposes Algorithm l6. ll resets only finitely-many times, e.g., 
after some iteration K A the sequence of problems enters the preasymptotic regime and 
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iterations on subsequent refinements are terminated by either the first or second condition 
of Criteria [6T2l 

Assumption 6.5. There exists a refinement K A , for which either condition (1) or condi¬ 
tion (2) of Criteria \6.2\ is used to terminate the iterations on refinement k, for all k > K A . 

Regardless of whether the exact problem within the target class is well posed, the 
approximate problems may not be. The coarse mesh problems are often ill-posed, may 
feature indefinite and ill-conditioned Jacobians, and there is no guarantee that even an 
exact solution to a given approximate problem will land in any basin of attraction for the 
approximate problem on the subsequent refinement. 

An analysis of the problem class where Algorithm 16.11 can be proven to converge, 
i.e., where Assumption 16.51 can be shown to hold, is beyond the scope of the current 
paper. The numerical examples in the next section demonstrate the performance of the 
algorithm on a range of problems with steep internal and boundary layers, as well as 
highly variable solution dependent frequencies in the diffusion coefficient. The following 
Theorem [6T6] summarizes the theoretical convergence properties of Algorithm l6.1 1 within 
the problem class characterized by Assumption l6.51 

Theorem 6 . 6 . Assume the hypothesis of Corollary 14.71 and Corollary 15.61 and let As- 
sumption 16.51 hold for k > Ka • Then the iterations on refinements k > Kb for Kb < 
K \ + Na + Nb with Na and Nb given by (16.41) and (16.51) . attain convergence up to 
tolerance £con with 5 = 1. Further assuming the hypotheses of Theorem 15.61 there is a 
refinement Kc for which 7 *. = 1 for k > Kc- 

The convergence of 7 to one is not necessary for the convergence of the iterations, 
but it is necessary for its efficiency, and the asymptotically quadratic convergence of the 
resulting iterations as discussed in [121 118112211211 . 

Proof The convergence of the sequence of residuals ||r fc || —y 0 within tolerance c CO n 
follows directly from Criteria 15711 Assuming ||rxJ| > £con and 7 < 7mono from 
Assumption 14.61 direct calculation shows after at most N A refinements, the iterations 
achieve convergence of the residual, with N A given by 

^ = r i°f wikid 1. (6 . 4) 

log(l — I/27MON0) 

In particular, ||ry | < l C on for k > K B — K A + N A - From Corollary 15.61 and Assump¬ 
tion [531 it also holds if 5k b < 1, the sequence 5/,. k > K B increases monotonically to 
unity, and for q 1 < q < 1 it follows that 5 k = 1 for k > I\ c = Ii B + N B , with 

Nb — Rog Sk b / log (fl • (6-5) 

This establishes the first result, the convergence of the sequence of residuals to zero for 
the problem defined with consistent data. 

The second part of the argument demonstrates the decrease of 7 to 1 follows from 
the local convergence Theorem 13.61 Assuming the residual is small enough, then the 
rate ||r n+1 ||/||r n || asymptotically approaches 1 — I/ 7 , and in particular is within e T 
tolerance of the predicted rate, meaning Criteria 14.21 hold, and 7 is updated by Defi¬ 
nition 14.31 and, by Corollary 14.71 converges to one after at most an additional N c = 
[— log(y M oNo)/ log q] updates. The subtle point is the convergence tolerance £con must 
be small enough for the asymptotic convergence rate of Theorem 13.61 to be assured. □ 

Remark 6.7. It is not just an artifact of the analysis that 7 —> 1 cannot be assured un¬ 
til after the iterations already achieve convergence to tolerance on every iteration, and 
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that the tolerance Econ must be sufficiently small. The phenomenon of 7 stalling above 
one and failing to update for several iterations after \ \ r k \ \ < e CO n has been observed, in 
particular in problems similar to Example 17.21 in the next section, where a diffusion co¬ 
efficient with a high frequency dependence on the solution fails to resolve until the mesh 
in the highest-frequency region is sufficiently fine. It is more commonly observed that the 
iterations do converge at the predicted rate when the residual is far from convergence, 
but generally 7 —> 1 within a few iterations on either side of the residual converging to 
tolerance. 


1. Numerical Examples 

The numerical experiments were performed with a Python implementation of the FEn- 
iCS library lfT9l . and additional Numpy IflTl computations. For each of the examples 
shown, Algorithm 16. II is demonstrated solving for both a known and an unknown solu¬ 
tion. The unknown solution examples demonstrate the effectiveness of the method for 
cases where the variable dependent and solution dependent layers may not match up. For 
the source functions designed to produce a known solution, the H l and L 2 errors denoted 
by | u — Uk |i and \u - ufa, respectively, are shown along with the error estimator; oth¬ 
erwise, the error estimator is shown to decay at the expected rate of n ” 1 / 2 for degrees of 
freedom n, after the resolution of the layers. 

The examples are all run with reduction factor q 1 = 0.9 and a 0 = 0.9. The Dorfler 
parameter of (13.31) is set to 6 = 0.2, at the tolerance Ccon = 10”'. Example 17.21 with 
the unknown source (17.81) uses rate tolerance £t = 0.01, and all other examples use 
Et = 0.005. The higher rate tolerance is used to prevent the update of 7 from stalling 
in cases where a high mesh resolution is necessary to accurately capture the Jacobian. 
Perturbing these parameter settings may change the number of refinements or the number 
of resets taken for the algorithm to arrive at the asymptotic regime, but not the end result. 
The initial and reset values of u, 7 , and 5 are set as in Remark 16.41 in accordance with 
Assumption 14.21 but not necessarily Assumption 14.61 Nonetheless, monotonicity of the 
parameters is generally observed on iterations when the algorithm does not reset. Each 
example is started on a uniform mesh of 144 triangular elements. 

Example 7.1 (Steep internal layer). Consider the quasilinear diffusion equation on Q = 
[0,1] x [0,1]. 

F(u,x ) := — div(K(w)Vw) — f(x,y) = 0 in f), u = 0 on dCl. (7.1) 

The solution dependent diffusion coefficient is given by 

ffis ) = k -1 - 7 -— with a = 0.5, k — 1, and e = 6 x 10” 5 . (7.2) 

£ + (s — a) z ' 

The following source functions are considered 

f(x, y) chosen so u(x, y) = sin (the) sin( 7 n/), (7.3) 

f(x, y) = 10 5 ■ (0.5 -\x- 0.51 ) 2 (0.5 - \y - 0.5|). (7.4) 

The diffusion coefficient k(s) is bounded away from zero with tf s), k'(s) and n"(s) 
all bounded assuring uniqueness of the solution as in @1. Milder versions of this model 
problem with source (17.31) are considered by the author in ll22l 21]. The inexact itera¬ 
tion (13.71) allows the solver to stabilize on a coarser mesh and with steeper gradients in 
the diffusion term than in previous versions of the method. Figure [Hon the left shows an 
adaptive mesh in the preasymptotic regime for Example 17.11 with source (17.31) showing 
local refinement of the internal layer with the mesh remaining coarse away from the layer. 
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FIGURE 1 . Adpative meshes from Example 17.11 Left: mesh after 20 iterations 
with 4140 elements for source (17.31) . Right: mesh after 25 iterations with 2208 
elements for source (17.41) . 


The mesh on the right shows the boundary and internal layer resolution for Example l7.ll 
with source (17.41) . 

Figure [2] shows the II ] error, the L 2 error and the error estimator running from the 
course mesh regime through the preasymptotic and into asymptotic regime. The transi¬ 
tion from the coarse mesh to preasymptotic is captured by the increase in error within the 
first twelve iterations, and through the nonmonotonicity of parameters 7 and 5, shown on 
the right. The decreases in 5, and the corresponding increases in 7 during this phase are 
due to reset of the algorithm. The transition from the preasymptotic to the asymptotic 
regime is displayed as a sharp decrease in L> error and the estimator 77 as 5 approaches 
one and the iteration solves for a right hand side consistent with the problem data. This 
sharp drop can be smoothed over several iterations by allowing a few preliminary refine¬ 
ments based only on the variable dependent data or otherwise stalling the initial increase 
of 5, but these techniques were not used in the current work in order to better illustrate 
the range of possible behaviors of the inexact method. In the plot of the parameters 7 
and 5, Figure [2] also illustrates how the initial relationship 7 = 1/5 is approximately 
maintained when the iterations maintain stability through the mesh refinements. 

Figure [3] shows the behavior of the error estimator next to a plot of the parameters 
7 and 5 for the diffusion problem of Example 17.11 with a nonsmooth source (17.41) . The 
initial coarse mesh problem is a bad representation of the exact problem resulting in 
early resets of the algorithm, after which enough of the data is captured to maintain 
stability through the mesh refinements. In contrast to the same diffusion operator with 
source (17.31) . here the error estimator grows throughout the preasymptotic regime as 5 
increases to one, then decreases at the predicted rate of with respect to degrees of 
freedom n. 

Example 7.2 (Oscillatory diffusion with several layers). Consider the quasilinear diffu¬ 
sion equation on Q = [0,1] x [0,1]. 

F(u, x ) := — div(/c(u) Vti) — f(x, y) — 0 in O, u = 0 on di1. (7.5) 

The solution dependent diffusion is given by 

k(s) = k + sin(s/e) + arctan(s/£), with e = 6 x 1CT 3 


(7.6) 
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FIGURE 2. Left: H 1 error, L 2 error, and error estimator (above) against n '/ 2 
where n is the number of degrees of freedom. Right: parameters 7 and 6 for 


Example 17. II with source (I7.3I ). 




FIGURE 3. Left: Error estimator 7 against n 122 where n is the number of de¬ 
grees of freedom. Right: parameters 7 and 5 for Example 17. II with source (17.41) . 


The following source functions are considered, with the k from (17.61) set to k = 2 + tt/2 
with source (17.71) and k = 10 with source (17.81) . 

f(x, y ) chosen so u(x, y ) = 1.85 4 • x(x — 1 )y(y — 1) — 2 8 • x 9//4 (a; — 1 ) 2 y 2 (y — l) 2 , 

(7.7) 

f(x, y) = (1 -x)(l -y)(e^ - l)(e 8 »* - 1). (7.8) 

By the same reasoning as Example 17.11 this problem has a unique solution. Exam¬ 
ple 17.21 demonstrates the inexact method on a problem where in addition to steep gradi¬ 
ents, the diffusion coefficient features oscillations proportional in frequency to the mag¬ 
nitude of solution u. Similar to the first example the error decreases throughout all iter¬ 
ations with minor fluctuations for the source with the known solution (17.71) . In this case 
for the source corresponding to an unknown solution (17.81) . the estimator first decreases 
then stabilizes and increases as S increases to one, and then decreases at a steady rate. 

For both source functions in Example 17.21 the decrease in 7 stalls for several iterations 
in the preasymptotic regime. This shows the iterations are not converging at the expected 
rate for several refinements, presumably because the solution-dependent frequencies of 
the diffusion coefficient k(u), are not be sufficiently well expressed on the mesh as the 
solution u first attains a threshold frequency. 
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FIGURE 4. Left: H 1 error, L 2 error and error estimator (above) against n -1 / 2 
where n is the number of degrees of freedom. Right: parameters 7 and 6 for 
Example 17. 2 1 with source (I7.7I ). 




FIGURE 5. Left: Error estimator against n l,/2 where n is the number of degrees 
of freedom. Right: parameters 7 and 5 for Example [L2] with source (17.81) . 


Example 7.3 (Gradient-dependent diffusion coefficient). Consider the quasilinear diffu¬ 
sion equation on = [0,1] x [0,1], 

F(u,x ) := — div(«(| Vu\ 2 )Vu) — f(x,y) = 0 in 0, u = 0 on dCl. (7.9) 

The solution dependent diffusion coefficient is given by 

n(t 2 ) = k + arctan((f 2 — a)/e), with a = n, k = n, and e = 2 x 10~ 2 . (7.10) 

The following source functions are considered 

f(x,y ) chosen so u(x,y) = sin( 7 r.x') sin( 7 rt/), (7.11) 

f(x, y)= 2 x 10 3 • (0.5 -\x- 0.5 |) 2 (0.5 - | y- 0.5 |) 2 . (7.12) 

As in fl4l . monotonicity and Lipschitz conditions are established by observing there 
are positive constants c and C with c < n(t 2 ) + 2 t 2 n'{t 2 ) < C, assuring a unique solu¬ 
tion to the exact problem. This is not known for the approximate problems, particularly 
on the coarse meshes. This example with source (17.1 II) illustrates the phenomenon of 
multiple resets as thin layers in both the source f(x) and the Jacobian are uncovered. 
Figure [ 6 ] corresponding to Example 17.31 with source (17.1 II) shows several jumps up in the 
error corresponding to the resets, a steep decrease as 5 —» 1 and the discrete problem 
regains consistency, then steady decrease in the error at the predicted rate. The plot on 
the right shows 7 stalling for several refinements after the last reset until the approximate 
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FIGURE 6. Left: H 1 error, L 2 error and error estimator (above) against n -1 / 2 
where n is the number of degrees of freedom. Right: parameters 7 and 6 for 
Example 17. 3 1 with source (17. ill ). 




FIGURE 7. Left: adpative mesh from Example 17.31 with source (17.121) after 12 
refinements with 1914 elements. Right: parameters 7 and 5. 

Jacobian is sufficiently well resolved and Criteria l4~2l for updating 7 are met before Cri¬ 
teria [5J] determines the iterations should stop. Once necessary resolution is achieved, 7 
reduces quickly to one and the residual is solved to tolerance within a few iterations on 
subsequent refinements. 

The decrease in the estimator is relatively steady in comparison to the H 1 error, and 
does not indicate that preliminary to the resets the correction steps w n fail to update 
the solution u n to sufficiently decrease the residual ||r n ||. In this case the update steps 
w n display a loss of stability due to a qualitative change in the problem data from the 
latest mesh refinement, albeit at a magnitude small in comparison to u n , and the reset is 
performed before the solution is dominated by unstable fluctuations. 

For source (17.121) corresponding to the unknown solution, the estimator decreases sim¬ 
ilarly to the one for source (17.1 II) . The adaptive mesh for the twelfth refinement, just as 
7 — » 1 is shown on the left of Figure [7] next to a plot of parameters 7 and 5 . The mesh 
shows the most dense refinement in the center, corresponding to the peak of the source, 
with additional refinement along the discontinuities of its gradient. Additional dense 
refinement is seen developing along the steep gradients of k(\ \7u\ 2 ). 

8 . Conclusion 

An updated adaptive method for the numerical solution of quasilinear diffusion prob¬ 
lems with steep internal layers is presented, building on the previous work by the author 
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in lf22ll2T1l . The goal of this method is to start the computation on a coarse mesh where 
the solution and variable dependent problem coefficients are not resolved, and to par¬ 
tially solve the possibly ill-posed approximate coarse mesh problems up to sufficient 
tolerance to refine the mesh and eventually achieve a stable and efficient computation 
of the solution to an accurate representation of the problem. The current results include 
explicit definitions for the added numerical dissipation parameter 7 used to stabilize the 
iterations on the coarse mesh, and the introduction of an inexact method to further sta¬ 
bilize problems where the residual is dominated by a variable dependent source. The 
scaling parameter S of the inexact method is shown to reach one which demonstrates 
consistency of the method, and the convergence of the residual of the discrete nonlin¬ 
ear problem is established. Under convergence of the residual, the parameter 7 is also 
shown to reach one, establishing efficiency of the Newton-like iterations in the asymp¬ 
totic regime. Future work by the author will include adapting the current method for 
quasilinear convection problems and nonhomogeneous boundary conditions, as well as 
investigation of the stabilization properties of inexact linear solves. 
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